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1. Motivation and objectives 

The most basic result in a study of decaying isotropic turbulence is the evolution 
of the kinetic energy as a function of time. By postulating a self-similar decay of 
the energy spectrum based on an exact invariant Bo of the flow, Saffman ( 1967a, b) 
determined the high Reynolds number decay law 

(u 2 ) oc B$t-i (1) 

where Bq is the leading coefficient of the energy spectrum near k = 0 

E(k) ~ 2nBok 2 k — ► 0. (2) 

Saffman ’s determination of the high Reynolds number decay exponent was based 
on earlier work by Kolmogorov (1941) in which it was assumed that a self-similar 
decay of the spectrum could be based on the invariance of the Loitsianski integral 
B 2 (Loitsianski, 1939), yielding the decay law 

(u 2 ) oc Bft^^ (3) 

where now 

E(k) ~ 27tB 2 k 4 k — ► 0. (4) 

However, it was later shown (Proudman & Reid, 1954; Batchelor & Proudman, 
1956) that B 2 was in fact not invariant and depended on time during the turbulence 
decay. 

However, one may still postulate an exact self-similar decay of the energy spec- 
trum at large-scales (Lesieur, 1990). If it is assumed that 

B 2 (t) = /3f>, (5) 

then (3) still holds but with B 2 (t) given by (5). When 7 is positive, as is indicated 
by numerical simulations and quasi-normal closure models, this results in a less 
rapid decay of the energy as £~ 10 / 7 + 2 7/ 7 

We have performed large-eddy simulations of decaying isotropic turbulence ( Chas- 
nov, 1994) to test the prediction of self-similar decay of the energy spectrum and 
to compute the decay exponents of the kinetic energy. In general, good agreement 
between the simulation results and the assumption of self-similarity were obtained. 
However, the statistics of the simulations were insufficient to compute the value of 
7 which corrects the decay exponent when the spectrum follows a k 4 wavenumber 
behavior near k = 0. To obtain good statistics, it was found necessary to average 
over a large ensemble of turbulent flows. We report on this work here as well as in 
a recent Physics of Fluids A letter (Chasnov, 1993). 
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2. Computation of the Loitsianski integral 

The coefficient B 2 above, the so-called Loitsianski integral, can be written as an 
integral over the infinite flow volume as 

52 = _ irrJ im /(“.(*)“<(* + r))r 2 dr. (6) 

487T’ V-*ooJy 

To compute B 2 by numerical simulation, we assume that the velocity field is 
periodic in three directions with periodicity length L = 2ir. The velocity field may 
then be expanded in a Fourier series as 

u(x) = ^ u(k) exp(tk • x), (7) 

k 

where the components of k in the sum span the set of integers. A good approxi- 
mation to homogeneous turbulence is thus obtained when the integral scale of the 
turbulence is much less than n. Treating the average in (6) as a volume average, 
and substituting the Fourier expansion (7) into (6), we obtain after one integration 
over the volume 

i?3 = ~ 4^7 ^ Ui(k)tij(— k) J exp(ik • r)r 2 dr. (8) 

k " 

The remaining volume integral in (8) may be evaluated analytically, and making 
use of Uj(— k) = u,(k)*, where * denotes the complex conjugate, and Ui(0,0,0) = 0 
we obtain 


b 2 = £ { S [|u(M,0)| 2 + |u(0,M)| 2 + |u(0,0,fc)| 2 ] . (9) 

There are two main difficulties in the direct use of (9) to compute B 2 in a numer- 
ical simulation. Firstly, the correlation (uj(x)ui(x + r)) in (6) decreases in general 
as 0(r“" 5 ) in homogeneous turbulence (Batchelor Sz Proudman, 1956) - although 
it decreases faster as o(r“ 6 ) in an isotropic turbulence - so that the integral scale 
of the turbulence must be small enough for the integral in (6) to converge within 
the computational domain. Secondly, as the value of r in (6) becomes comparable 
to 7 r, the replacement of the ensemble average in (6) by a volume average becomes 
inaccurate because of a lack of sample of the largest computed scales. Explicit 
computation has shown that direct use of (9) to compute B 2 in a single realization 
of a turbulent flow is highly inaccurate. We are thus led to average B 2 over an 
ensemble of such flows. This is equivalent to treating the original average in (6) as 
a combination of a volume and ensemble average. 

In this research brief, we report on a computation of B 2 (t ) accomplished by 
performing 1024 independent simulations of resolution 64 3 . The size of this en- 
semble is sufficient to compute B 2 (t) to a statistical accuracy better than 5% over 
the entire time-evolution considered. The computations are performed on an Intel 
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iPSC/860 hypercube parallel machine containing 128 processors. The machine had 
eight megabytes RAM per processor which allowed 64 realizations to be performed 
in parallel with each independent realization computed on 2 processors. Commu- 
nication between processors computing different realizations is minimal so that the 
simulation of an ensemble of turbulent flows makes very efficient use of parallel 
computer architectures. Sixteen independent runs — each of 800 total time-steps 
— were performed. With each time-step taking approximately 10.6 seconds of cpu 
time, a total of 38 hours of dedicated machine use was required. 

Our main goal in computing i? 2 (<) is to determine its long time, high Reynolds 
number behavior. Under the constraints imposed by 64 3 resolution simulations, this 
necessitates the use of a large-eddy simulation with the initial peak of the energy 
spectrum placed at as large a value of k magnitude as possible (Chasnov, 1994). 
Here, the initial energy spectrum is taken to be 

E(k,0) = 27r£ 2 (0)fc 4 exp \-2(kjk p ) 2 , (7) 

with k p = 25 and S 2 ( 0) = 6.934 x 10~ 8 , so that (u 2 ) — 1. As we have done previ- 
ously, an eddy-viscosity subgrid scale model (Kraichnan, 1976; Chollet k Lesieur, 
1981) is used to model the unresolved small-scale turbulence. Although the inclu- 
sion of a stochastic backscatter term in the subgrid model (Chasnov, 1991) can 
directly affect the time- variation of B 2 , this effect is negligible at the later times of 
the turbulence evolution of interest to us here. 

The finite resolution of the simulation results in a spherical truncation of the 
Fourier series in (7) at k m , the maximum wavenumber of the simulation, so that 
the sum to co in (9) is replaced by a sum to k m j\J 3. At small times when the 
peak of the energy spectrum is near k m , this sharp cutoff results in errors in the 
computed value of f? 2 . We have shown that these errors can be easily removed by 
applying an additional Gaussian filter of the form exp[(-Jt/Jfc/) 2 ] with k f = 12 to 
u(k) before computing (9). At the later evolution times of interest to us here, the 
effect of this additional filter is negligible. 

The results obtained from the simulations are shown in figures 1-3. In figure 1, we 
plot the time-evolution of the ensemble- averaged energy spectrum obtained from the 
large-eddy simulations by summing the contributions of |u(k)| 2 into wavenumber 
shells of thickness Ak = 1 in the usual way, i.e., 

2nk 2 

£(M) = -£— 

where S* is the number of Fourier modes in each wavenumber shell and k = 
1.5, 2.5, 3.5, 29.5. A good approximation to the homogeneous turbulence en- 
ergy spectrum is thus obtained at high wavenumbers, while the approximation is 
less accurate at low wavenumbers. Nevertheless, the increase in time of the low 
wavenumber k 4 coefficient is clearly evident from the plot. 
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FIGURE 1. Time-evolution of the energy spectrum. 



*/t(0) 

FIGURE 2. Time-evolution of the Loitsianski integral. 

The coefficient B 2 (t)/B2(0) versus time, in units of the initial large-eddy turnover 
time r(0) where r(0) = 1.38/(fcpJ32(0)) 1 ^ 2 , is plotted in figure 2. The points rep- 
resent the statistical mean of the ensemble while the pluses represent one standard 
deviation from the mean. The standard deviation of the distribution of B 2 itself, 
which we have shown from the simulation data to be approximately Gaussian, varies 
somewhat over the course of the simulation but at the latest time plotted is 80% 
of the mean. With 1024 realizations, the statistical uncertainty of the mean at the 
latest time is 2.5%. 
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FIGURE 3. Time-evolution of the logarithmic derivative of the Loitsianski integral. 

In figure 3, we plot the logarithmic derivative of B 2 with respect to time in order 
to determine the validity of (5) and to compute a value of 7 from the simulation. 
In agreement with the EDQNM model, we find that ^(0 follows an approximate 
power-law at large times. From figure 3, we estimate the power law exponent to be 
7 ~ 0.25, with a statistical uncertainty of 6 % at the latest time. The straight line 
drawn on the log-log plot of figure 2 represents this result. The value of 7 we obtain 
from the simulation is about 50% larger than that estimated previously (Lesieur &: 
Schertzer, 1978; Lesieur, 1990). Using our computed value for 7 , the Kolmogorov 
decay exponent becomes —1.36 instead of —1.43, a difference of 5%. 

The statistical uncertainty of our asymptotic result for 7 can be reduced further 
by computing additional realizations. However, there may be other errors in our 
result associated with the deviation of “periodic turbulence” from homogeneous 
turbulence at the latest times of evolution, as well as the expected slow approach 
of the turbulence to asymptotics (Chasnov, 1994). The evident trend of figure 3 
is towards a somewhat smaller asymptotic value for 7 than we have estimated. It 
would be of interest to repeat the present computation at higher resolution with a 
larger ensemble after parallel machines have become substantially more powerful. 

We also note here another approach to the current computation. Rather than 
simulate 1024 64 3 turbulent fields, we could have simulated 16 256 3 fields with 
slightly more computer time due to the need for inter-processor communication. To 
obtain similar statistics between these two simulations, we would have to increase 
the initial peak of the energy spectrum k p by a factor of four and truncate the 
volume integration in ( 8 ) to 1/64 the volume of the entire periodic domain. It is 
unclear which simulation would result in a more accurate computation of l? 2 (t), 
but we chose the former mainly to illustrate the efficiency of performing realization 
averages of turbulent flows on parallel machines. 
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3. Conclusions 

This work has demonstrated the capability of numerical simulations to compute 
large-scale statistics of turbulent flows by means of an ensemble average over a large 
number of independent realizations of the flow. Such a technique is “embarrassingly 
parallel” and is ideally suited for the new parallel computer architectures. This tech- 
nique may also be applicable to turbulence simulation on virtual parallel machines 
in which many powerful workstations are connected together over a local network. 
If the memory of each workstation is sufficiently large so that each realization can 
be performed independently on each workstation, then the only communication 
required between workstations is to perform the ensemble average. 
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